function vals_ = ddPow6(fits, data)
% function vals_ = ddPow6(fits, data)
%
% Computes TIME-based DDleak function.
%   Assumes DD to a fixed TIME (given as col 2 of data).
%   Thus, pct correct is fraction of gaussian above
%       0 at TIME.
%   Drift rate depends on coherence as a power law:
%       u(coh) = A*coh^M
%   Bias is fixed offset plus noise
%       bias = normrnd(B, Bsd)
%
%   Accumulator leaks with "leakiness constant" L. Thus
%   mean  u(t) = bias_mean + u/L*(1-exp(-Lt))
%   SD   sd(t) = sqrt(bias_var + FF *mean/L*(1-exp(-2*Lt))))
%       where
%           FF is fano factor, var/mean for S1 and S2
%
%   at values in "data":
%       data(1)   ... coh (%)
%       data(2)   ... time (s)
%       data(3)   ... dot dir (-1/1)
%
%   given parameters in "fits":
%       fits(1) ... A      (coh scale)
%       fits(2) ... M      (coh exponent)
%       fits(3) ... N      (time exponent)
%       fits(4) ... B      (bias)
%       fits(5) ... Bvar   (var of bias)
%       fits(6) ... lambda ("lapse")

% return initial values (init, min, max)
if nargin < 1
    
    vals_ = [ ...
        30    0.1   100;   ...
        1     0.001  10;   ...
        1     0.001  10;   ...
        0   -10      10;   ...
        0     0   10000; ...
        0     0       0.45];

else
    
    R0 = 10;
    S1 = (R0 + fits(1).*data(:,1).^fits(2)).*data(:,2).^fits(3);
    S2 = R0.*data(:,2).^fits(3);
    
    mu = data(:,3).*fits(4) + (S1 - S2);
    sd = sqrt(fits(5) + 0.3.*(data(:,3).*fits(4)+S1+S2));
    sd(sd<0) = 0;
    
    % compute pct correct from the dvar
    p = 1 - normcdf(0,mu,sd);
    
    % Avoid NaN's in normcdf, which happens when s=0 (i.e., time=0)
    p(isnan(p)) = 0.5;     % otherwise is NaN
    TINY        = 0.0001;
    p(p<TINY)   = TINY;    % lower bound
    p(p>1-TINY) = 1-TINY;  % upper bound

    % abbott's law
    % for dealing with performance that does not asymptote at 1
    % C = performance at chance
    % D = asymptote as stim -> inf
    % Abbott's law: p* = C + (1 - C - D)P
    % where P is percent correct on [0 ... 1]
    % But remember that the p we just computed is [C ... 1], C = 0.5;
    % i.e., p = C + (1 - C)P
    vals_ = 0.5 + (0.5 - fits(6)).*(2.*p-1);
end
